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Abstract- This paper deals with the comparison of three types of 
sub-optimal control law for the stable levitation of a turbojet 
shaft, sustained by two radial active magnetic bearings (AMBs). 
Shaft is considered rigid for the procedure simplification. The 
utilized approach leads to development of different sub-optimal 
control laws to use in speed-varying simulations in the angular 
speed of the shaft. The first control matrix is obtained by explicit 
relationships of the parameters of the control law vs. speed, 
obtained using a curve-fitting procedure neglecting the speed- 
varying elements out of the main diagonal of each single block 
constituting the entire control matrix. The second control law is 
obtained by control matrix set to zero angular speed during 
speed-varying simulations. The third control matrix is obtained 
by the second neglecting the elements out of main diagonal. Time 
response to the step input shows difference in performance of the 
AMBs when varying the control law. All control laws are tested 
by means of simulation. The study allows showing like a 
polynomial representation of the control matrix, developed to the 
variation of angular speed, is able to maintain stable the 
magnetic levitation. 

Keywords- Optimal control; Sub Optimal control; Active 
magnetic bearing; Turbojet shaft; Closed loop 

I. INTRODUCTION 
The system considered in this paper consists of a floating 
mechanical rotor and three active magnetic bearings (AMBs), 
two radial and one axial, which provide suitable controlled 
dynamic forces, thus allowing the suspended object to display 
a stable motion. 

Due to the contact-less operation capability, the AMB 
system has many promising advantages for high-speed, high- 
temperature, and clean environment applications. Moreover, 
adjustable stiffness capabilities and damping characteristics 
also make the AMB suitable for elimination of vibrations. 
Although the system is complex and considered an advanced 
topic in terms of its structural and control design, the 
advantages that it offers compensate for the design 
complexity. 

Application areas of magnetic bearings are still steadily 
expanding. Some AMB applications that have received many 
attentions [22] by many research groups are the flywheel 
energy and storage device [l]-[2], the turbo molecular pump 
[3], compressors [4], the Left Ventricle Assist Device (LVAD) 
[5] and the artificial heart [6], [7], [8]. 

Magnetic bearing systems are designed in a few 
configurations to satisfy various specifications that arise from 
different applications. The rotor can be oriented in horizontal 
and vertical positions; in both positions, the effect of gravity 
is decoupled from the system dynamic forces and, 



consequently, the vertical displacement is controlled 
separately from the other set of magnetic coils controlling 
horizontal position [3], [9]. 

The magnetic bearing system consists of a rotor in the 
horizontal orientation. However, it is, more widely used. 
Moreover, in some applications, such as the artificial heart, 
the rotor is oriented in both the horizontal and vertical 
positions, due to the nature of operations depending on the 
position and movement of the host [7], [8]. 

The design of magnetic bearings is carried out utilizing 
suitable methodologies. Often, a linear model of these 
actuators is used to approximate the nonlinear relationship 
between force, current, and air gap length, thus obtaining a 
model suitable for the design of linear controllers. 

A disadvantage of this approach is that the linear model is 
approximated around a single operating point, and the validity 
of this model decreases as the physical system is perturbed 
from this point. 

This approach has been previously considered for an 
AMB system in [10]. Other related approaches, such as input- 
output linearization [11] [12] have also been examined. In 
[13], experimental results from a one-degree-of-freedom 
AMB are given. In [14] a control law is derived and tested 
experimentally. In particular, it is shown that instability 
phenomena on rotor positions can exist due to coupling 
effects that result fromthe elements of the control matrix. The 
rotor dynamic modelling and experimentally extracted 
transfer functions are presented and analysed in the paper [21], 
[22]. 

However, if good feedback control is designed, the 
variations of the mechanical variables, such as the 
displacements of the shaft within the air gap between stator 
and rotor can be always kept very small and, consequently, 
the above linearization approach yields valid results. 

A stable electromagnetic suspension system may be 
obtained by means of state feedback optimal control. Because 
a turbojet shaft is subjected to angular speed variations during 
operations, a speed-varying control law is necessary. This 
control law should be obtained by solving, on-line, the matrix 
Riccati equation for each value of the angular speed, but this 
requires considerable computational efforts. 

The contribution of this paper compared to the previous is 
in the fact that the speed-bearing optimal control law is 
obtained by solving, off-line, the matrix Riccati equation for 
each value of the angular speed, thus obtaining a set of gain 
matrices. In succession, the analytical expressions of all the 
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elements of these gain matrices vs. the shaft speed are 
determined using a curve-fitting procedure. Several 
simulation tests carried out on a closed-loop system consisting 
of the rotating shaft and the AMBs show that this optimal 
control law, figured in polynomial matrix form, can support 
the rotor at every value of the angular speed in presence of 
both static and harmonic loads. No experimental parts are 
introduced with the present paper, because the experimental 
tests were executed in the paper [24], with the purpose to 
validate the theory presented in this paper. 

Hypothesis of rigid rotor is done due to the simplicity of 
the analysis; computational analysis of flexible rotor 
configuration has been left out for future developments in this 
field. 
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Fig. 1 Active Magnetic Suspension 




Fig. 2 Scheme of the Rotating Shaft 
Supported by Two Radial Active Magnetic Bearings 

II. MATHEMATICAL MODEL 

The basic scheme of an active magnetic suspension is 
given in Fig. 1, wherein a rotating shaft is maintained 
according to a well-defined configuration by means of three 
magnetic bearings, two radial and one axial. Rotor shaft 
parameters are reported in Table I. The control of the 
configuration is carried out modulating the current supply to 
the magnetic bearings. Figure 2 shows the shaft, the centre of 
mass, the forces generated by the radial bearings, reference 
axes x, y, z, the Euler's angles 3,<p,y and the reference 
sections A and B. 

Since the radial behavior is independent on the axial 
behavior, there is no reason to complicate the mathematical 
model to implement the control law that stabilizes the radial 
position of the shaft. 



Assuming that the shaft is a rigid body, the mathematical 
model of the whole system under analysis, with reference to 
the center of mass, is given by: 

Table I Rotor Shaft Parameters 



Symbol 


Quantity 


SI 


m 


Mass of rotor 


144 kg 


Jr 


Inertia along x-y 
direction 


26,87 kgm 2 


Ja 


Axial inertia 


1,65 kgm z 


£ 


Maximum value of the 

eccentricity of the mass 

center 


4 10*m 


X 


Principal inertia axis 
deviation 


6 10* 


n 


Contingency coefficient 


13 



M i5 + O.G^ = Bi(i,q ) + f +f(Qt) 



dt z 
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Considering the coordinate transformation: 

q b =kA y bA x b B y bB f =K^ 



(2) 



the following model involving the displacements of sections 
A and B is obtained: 



M b q b + QG b q b = f b + nB q \ + B^f (Qt) 
where the components of f . are given by: 



(3) 
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IxB = K 
fyB = K B 
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(3a) 



where: 



k A = 1/4 jU NlA A cos a Jc B = 1/4 ju Q N 2 B A B cos a 

n is a contingency coefficient and 2a is the angle between 
two successive magnetic poles (a=22°,5 in this work); g A and 
goB are the air gape of the two bearings A and B respectively 
(their values are reported in table II); ii are the total currents in 
the coils of the electromagnet (they are four for every bearing). 
In Aeronautics the contingency coefficient "n" indicates how 
many times it must multiply the acceleration of gravity to 
obtain the acceleration of the airplane in the curvilinear path. 
In this work, as the rotor is part of the airplane, then it is 
subject to a contingency coefficient of force acting on its 
center of mass. 

M b and G b are computed from M and G taking into 
account the change of variables: 



M„ 



Jr + lbB 2 m 
(ibA+lbB) 2 


-Jr + lbAlbBm 



(IbA+lbB) 2 




Ja 
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(IbA+lbB) 2 



-Jr + lbAlbBm 
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(IbA + lbB) 2 




(3b) 



(3c) 



Model (3) is nonlinear for the presence of magnetic forces 
f b . By means of linearization of these forces, the following 
model is obtained: 

M b q b (t) + QG b q b (t) = K s q b (t) + K,i A (t) + nBjjf g + Bjjf (Qt) (4) 
where: 

K s=< ii og(k xA ,k yA ,k xB ,k yB ) 



K, = diag(k ixA ,k iyA ,k ixB ,k iyB ) 



kxA = k yA = 2k A ^f- 
9oa 



k — k —Ik -^ 

K xB - K yB - ZK B 



3 

9ob 



k — k —7k 0A 

K ixA — K iyA — ZK A 



9oa 



k -k —7k -^ 

K ixB — K iyB — ZK B 



Sob 



Table II Parameters of the control law 



Symbol 


Quantity 


SI 


kxA - kyA 


Displacement gain of 

active magnetic 

bearing A 


17MNAn 


kxAB = kyB 


Displacement gain of 

active magnetic 

bearing B 


7,2 MNAn 


kixA - hyA 


Current gain of active 
magnetic bearing A 


5,6 kN/A 


kixB - hyB 


Current gain of active 
magnetic bearing B 


2,3kN/A 


J-maxA ~ J-rruxB 


Maximum value of 
current 


13A 


lOA - Job 


Bias current 


6 A 


Qoa- Qob 


Air gap 


2 10""m 



The static load on bearing B turns out to be smaller than 
the static load on the bearing A, because it is more distant by 
the center of mass. The proportionality constant of the 
magnetic force ka and kb are placed in the numerator, then the 
gains of position and current of A are larger than B. Of course, 
these gains affect the control law because the gain of position 
premultiplies the four position components of the state vector, 
influencing the values of the array A; instead the array of 
current gains premultiplies the vector of inputs (current 
control), representing the input-state transition array B. Both 
the arrays A and B are essential to formulate the stabilizing 
solution of Lurie-Riccati for the formulation of the optimal 
control law. 

As it is easy to verify, model (4) is unstable and 
consequently a suitable stabilizing control law has to be 
constructed. Putting: 



*(0 = [qj(0 Qb(t)] 1 



(4b) 

the state space form of (4) is given by: 

z(0 = Az(t) + Bi A (t) + B d [f g+ f(Qt)] 

q b (0 = Cz(0 
where: 



(5) 
(6) 



M-'K S -M-'G b Q 
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(6a) 



B, 



(4a) 



M b B 



, C = [I 0] 
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III. OPTIMAL CONTROL LAW 

In this paper we design a LQR control law of the shaft 
position. More precisely, the formulation of the control 

problem is: compute the control law i A (t) that minimizes 

the performance index: 

00 

J(«(0.iA(0)«„ = J[*(0 , 'Q«(0+^(0 r WA(0]* (7) 

o 
where z(t) = Qb(0 Qb(0 ' su bJ ect to the constraint 

(5)-(6); Q G R x and R e R 4x4 are, respectively, positive 
semi-definite and positive definite matrices obtained by mean 
Bryson's rule according to [23] the matrix Q and R are 
respectively: 



Q = diag < 



R = diag 



(maxx M ) (maxy M ) (maxx bB ) (maxy bB ) 

1 1 1 1 



1 (8) 

,0,0,0,(H 



(9) 



( max UJ (maxz^) (maxz cBx ) (maxz cBy ) 



where 



max x bA = max y bA = max x bB = max y bB = 0.001 m 



max z^ = max i cAy = max i cBx = max z cBy = 13 A 

Entries in the Q matrix, corresponding to the velocities 
when designing the LQR, are zero because the filling rule 
dictated by the references [14] is applied. 

The solution of the problem is given by [15], [16]: 

i A (t) = -K-z(t) (10) 

where K = R~ B S , S is the solution of the matrix 
Riccati's equation: 

A r S + SA-/?SBB r S = (11) 

The expression (11) refers to algebraic Riccati's equation. 
In this equation we can see the matrix of the dynamic system 
depending on angular speed that causes coupling effects about 
displacements in radial directions. If we have to build up a 
speed-variant control, a continuous computing of stabilizing 
solution of S matrix must be taken into account. This could be 
very burdensome from a computational perspective because 
the controller must solve in few seconds a system of 

n(n + l)/2 equations. 

This paper considers a different approach to develop a 
speed-variant optimal control the matrix K is obtained once 
the value of angular speed is known without solving that 
system of equation. This is obtained by a Matlab routine 
calculating each coefficient varying with the angular speed 
whose equation is carried out by a least squares curved fit as a 
second-order polynomial expression (it is executed clicking 
on the chart -Basic Fitting Tool and choosing the desired 
approximation and clicking on show equation). This leads to a 
control matrix K that is characterized by two block matrices 
(12), each is an algebraic matrix in the angular speed variable. 



k = [k k, 1 

L % % J 



where: 



K n =Q 2 A +QB.+C 

q b a b c 



K, =Q 2 D. +QE +F f 

q b d e f 

4x4 



(12) 



(12a) 



Where A a , B b , C c , D d , E e , F f E R 

Figs. 3 and 4 show the shapes of two of the 32 elements of 
matrix K, together with their analytical expressions. These 
were obtained in the Matlab environment by means of a 
curve-fitting procedure. The results of this procedure, when 
applied to all the entries of matrix K, show that the analytical 
approximation of these entries is at most a second-order 
polynomial. It follows that matrix K displays the structure 
reported in [17]. The resulting control law is given by: 

i A (t) = -[K qb K, b ][q b T (t) ql(t)J (13) 



K(l,l) = - 4.9e-006V - 3.7e-006*x + 7.6e+003 
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Fig. 3 Entry k 11 vs Speed. The Red Line Denotes the Analytical 
Approximation. 
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Fig. 4 Entry /cl5vs speed. The Red Line Shows the Analytical 
Approximation. 
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Fig. 5 Campbell Diagram with Speed-variant Optimal 

Control by Polynomial Approximation 
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Fig. 6 Campbell Diagram Produced by a Control 
with Constant Gains at Rad/s 

The result is characterized by the same dynamic behaviour 
with a minimal computational burden of the controller as 
showed by plotting the Campbell diagram of the system for 
speed-variant optimal control by solving the ARE and one 
with constant gains set to Orad/s shown in the Figures 5 and 6 
respectively, [18]. In them J a is the polar inertia moment of 

the shaft and J r is its transverse moment. A constant control 

matrix can be used, because each sensor is located along each 
axis as to capture the displacements and speed signal along 
that axis. The elements out of the main diagonal of each block 
matrix constituting the entire matrix K are neglected to 
analyse the dynamic behaviour of the shaft during its rotation. 

This technique was used in [17], with a lower value of 
rotor's mass and a difference of two order in polar moment of 
inertia. Coupling effects were measured in respect to the 
direction of gravity, according with the same direction of 
static acting load, along y-direction applied to the center of 
mass. The results show that in the presence of speed-variant 
optimal control, a displacement occurs along the x-direction. 
This displacement is caused by an increasing value of the 
control current due to the coupling effects of gyroscopic 
matrix as in the following expression: 



Wb]>4b]>tor 
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Fig. 7 Campbell Diagram Produced by a Control Matrix Set at Orad/s and 
Characterized by Two Diagonal Block Matrix 
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Fig. 8 Displacement of the Geometric Centre of Section A along y-direction 

This means that the generic control current signal is 
generated by linear combination between the elements of 

matrix and components of the state vector z(t) that causes 

displacements along the x-direction as fig. (7), (8) show. 

In Figure 5, the step represented by the Campbell diagram 
is the same in Figures 6 and 7. The three steps in these figures 
say that there is a sudden change of rotational frequencies for 
the modes of precession. The step is smaller in Figure 5 
because the optimal solution is obtained minimizing the 
values of angular velocity for which the sudden changes of 
the dynamic behavior may lead to instability in the rotating 
system Both figures 6 and 7, considering a variation of the 
angular velocity of the rotor, are generated by a non-variable 
control law with the speed, then both the control laws are not 
the optimal solution of the model, and cause the rotor to 
transfer more energy during the sudden changes of the 
dynamic behavior; they are dictated by step highlighted by the 
eigen frequency, then it are at higher angular velocities. 
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The displacements along the orthogonal direction to the 
static acting load are not influenced by gyroscopic effects as 
happens in the case when variable gains are used. Control 
current signals show the same behaviour depending on 
variable or constant gains as shown in the figure (9) and (10) 
from [17], [18], [19] and [20]. 
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Fig. 10 Waveform of Currents Flowing in the Upper and Lower 
Electromagnets of Bearing A 

However, a variable gain is introduced to show a different 
way to compute speed variant optimal control that is 
characterized by a similar dynamic behaviour to speed variant 
optimal control obtained by solving ARE (Algebraic Riccati 
Equation) and the produced control signal is obtained 
introducing only a value of angular speed rather than solving 
a complicated system of equations. 
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Fig. 11 Waveform of Currents Flowing in the Upper and Lower 
Electromagnets of Bearing B 



This model of speed variant optimal control is more 
accurately named "sub-optimal" because of its similar 
dynamics with respect to the solving method of Riccati's 
equation, whose Campbell diagram is shown in figure 11. 

IV. SIMULATION RESULTS 
The closed loop model consisting of the model (5)-(6) and 
the control law (11) is now tested by means of simulation. 

The simulations are carried out by the introduction of step 
input. It is introduced as exogenous excitation indentified by 
the gravity force, This acts at a certain time value in this paper, 
as shown in all figures at t equal to 0,2 s; the reason of the 
used step input is justified by the rapid movement of the 
aircraft during the fast rump up. In fact during the rapid 
climbing of aircraft, all the masses constituting the aircraft are 
subjected to a variation of acceleration. This is identified by n 
time the gravity force, where n is contingency coefficient, a 
typical parameter used in the structures design. 

Table I shows that transverse moment of inertia is greater 
than one order than the polar moment of inertia. This data are 
obtained using a rotor having almost the same mass of the 
rotor installed in a Rolls Royce engine model used in military 
aircraft. 

Figures 5 and 6 show the shape of imaginary parts of the 
eigenvalues of the closed loop dynamic matrix A — BK as 
it changes with rotor speed. 

In order to test the effects of the control law (13), the 
model (4) is implemented in Matlab-Simulink software 
putting: 

mscos(flt) 
Z (J r -J a )cos(nt) (14) 

m^sin(Qt) 
- Z {J r -J a )sm(Qt) 

Figs. 8 and 9 show the waveforms of the radial 
displacements of the geometric centres of sections A and B 
along the same direction of the gravity force (y-direction) 
multiplied by n for simulating the rapid pull-up of the aircraft 
in presence of the harmonic load (12), assuming 
Q = 7 Ot rad/s for t e [0, 2] and Q = 1 4 rad/s for 
t>2. 

Examination of Figs. 8 and 9 reveals that sinusoidal 
oscillations arise during the shaft rotation along both x and y 
axes. The maximum displacement of the section centres 
occurs along y-axis due to the combined effects of the gravity, 
amplified by a contingency coefficient to simulate a fast pull- 
up of the aircraft, and harmonic load. However, both 
oscillation amplitudes and displacements are very small - less 

than 8xl0" 4 m. 

The waveforms of the currents flowing in the coil 
windings of the active magnetic bearings are shown in Figs. 
10 and 11. In each of the two figures there are two traces. A 
trace refers to the upper electromagnet of a bearing, whereas 
the other trace refers to the lower electromagnet of the same 
bearing. Examination of Figures 10 and 11 shows that in both 
bearings of sections A and B, the maximum current is less 
than the design value, equal to 13 A. 



f(Qt) = Q 2 
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Considering [14], by setting the matrix K at zero, angular 
speed is built up to create stable levitation to develop the 
control of the shaft away from unstable phenomena. In this 
case matrix K has the structure: 
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Fig. 12 Waveform of Currents Flowing in the Upper and Lower 
Electromagnets of Bearing A 

(sub-optimal control at Orad/s) 
Figures 12 and 13 show the control current pathway. 
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Fig. 13 Waveform of Currents Flowing in the Upper and Lower 
Electromagnets of Bearing B 

(sub-optimal control at Orad/s) 

Now a simplification is considered: matrix (15) is set 
neglecting the elements outside the main diagonal of each 

block R 4x4 so that to have matrix K structure as (16): 



K = 



x10 



Xi 

o 






rv 22 
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^15 







^26 



















(16) 
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Fig. 14 Displacement of the Geometric Centre of Section A 
along y-direction 

Results are shown in Figures 14 and 15 for the 
displacements, and Figures 16 and 17 for the control current. 

Table III 
Displacement and control current comparison 

by matrixk in the forms (12), (15) and (16) 





K(12) 


K(15) 


K(16) 


ybAmax 


-8 HTm 


-55 1(Tm 


-55 1(Tm 


ybBrrux 


-12 1(Tiii 


-3 10 " 4 m 


-3 HTm 


IcAmax 


13A 


10,5 A 


10,5 A 


IcBmax 


15A 


9 A 


9 A 



Table III shows that matrix (12) doesn't follow the 
condition of Bray son's rule because displacements are greater 
than one millimeter. 
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Fig. 15 Displacement of the Geometric centre of Section B along y- 
direction 
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Fig. 16 Waveform of Currents Flowing in the Upper and Lower 
Electromagnets of Bearing A 
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Fig. 17 Waveform of Currents Flowing in the Upper and Lower 
Electromagnets of Bearing B 

It also happens for the control current; in fact, they are 
more thirteen ampere for the radial active magnetic bearing B. 
On the contrary, matrices (15) and (16) have similar 
behaviours in regard to displacements and control current. 
They are capable of maintaining displacements and control 
current signals around the values established by Brays on's 
rule. From this first comparison it seems that a speed varying 
matrix K given by relationship (16) produces displacements 
and control signals away fromthe design's values. Matrices in 
the form (15) and (16) produce almost the same effects in 
terms of displacements and control currents. Now consider the 
transfer function matrix to analyse the performances of the 
turbojet shaft when it is subjected to the different control 
matrices (12), (15), and (16). So that output is calculated as: 

y(s) = W(s)u(s) (18) 

where u(s) is the vector of the step input. 

Varying the angular speed, the step response of the entire 
system with control matrix (12) is showed in Figure 18. 
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Fig. 18 Step Response of Transfer Function Matrix 

Only the elements along the main diagonal of the transfer 
function matrix show a relevant response to the step input 
varying with the angular speed. Others are small more than 
two orders as shown in Figure 19. Based upon the results 
shown in Figure 19 the same elements are neglected in order 
to obtain the rise time and settling time from the elements 
along the main diagonal. 
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Table IV 
Rise time (T y ) and settling time (T ze o /o) ) of response to the step input 

VARYING CONTROL MATRIX (12) (15) AND (16) 





K(12) 


K(15) 


K(16) 




Ty[5] 


0,00648 


0,00663 


0,00669 


Active Magnetic 
Bearing A 


T Z (2%) [s] 


0,0107 


0,011 


0,011 


Ty[5] 


0,00507 


0,00499 


0,00524 


Active Magnetic 
Bearing B 


Tz(2%) [s] 


0,00844 


0,00825 


0,00883 



Table IV shows that the radial active magnetic bearing A, 
controlled using the control matrix (12), produces a lower rise 
time and lower settling time than one controlled by matrix (15) 
and (16). With regard to radial active magnetic bearing B 
using the control matrix (12), bearing B has a lower rise time 
than one obtained using the control matrix (16) but greater 
than one using (15). The settling time is lower than the one 
using the control matrix (16) but greater than one using (15). 

Of course the simulation results should be proven by 
experimental data, but at present there is no plan of such a 
study. 

V. CONCLUSIONS 
The paper describes a speed-varying optimal control law 
developed for active magnetic suspension bearings used for 
sustaining a rotating shaft. Simulation shows that the 
proposed control laws are able to maintain the shaft in a stable 
levitation in the presence of static and harmonic loads. The 
simulations show that contactless motion between stator and 
rotor is guaranteed during the increasing of angular speed in 
presence of acting loads. Moreover, all the involved variables 
stay within their maximum design values except when using 
the first control matrix. Response to the step input shows that 
the radial active magnetic bearing A has a faster response and 
dynamic precision using the first control matrix rather than 
using the last two. The same can not be said for the bearing B, 
because it exhibits good performance in terms of rise time and 
settling at different control law, the settling is still faster than 
that of bearing A (table IV). This is because B is more distant 
from the center of mass and is therefore subject to a smaller 
static load, with a different dynamic behavior as a function of 
the control law. 

However, the radial active magnetic bearing B shows 
good performances depending on the control matrix used. The 



future work should include a) the comparison of this control 
law with other control laws in the topic literature; b) make a 
comparative study between symmetrical bearing positions and 
asymmetrical bearing positions; and c) the instability analysis 
depending to the control law. 
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